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Abstract 

t — ■ We consider the problem of nonlinear filtering of one-dimensional diffusions from noisy 

measurements. The filter is said to lose lock if the estimation error exits a prescribed 
| region. In the case of phase estimation this region is one period of the phase measure- 

5_! ■ ment function, e.g., [— 7r, it]. We show that in the limit of small noise the causal filter 

that maximizes the mean time to loose lock is Bellman's minimum noise energy filter. 

1 Introduction 

U" 

q | Optimal filtering theory defines different optimality criteria, such as minimizing the condi- 

tional mean square estimation error (MSEE), given the measurements PQ, maximizing the a 
posteriori probability (MAP) density function (pdf) of the signal, given the measurements 
[2], Bellman's criterion of minimum noise energy (MNE) [3], [I], [5], and more. In problems 
of phase estimation, that lead to loss of lock and cycle slips, an important optimality cri- 
terion is maximizing the mean time to lose lock (MTLL) or to exit a given region, which 
^ ; is also a well known control problem [6J, [TJ, [8], [9]. Approximation methods for finding 

the various optimal niters have been devised for problems with small noise, including large 
i/-) ■ deviations and WKB solutions of Zakai's equation, the extended Kalman filter (EKF) |10j . 
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[H][T2], [13] and others. The EKF and WKB approximations produce explicit suboptimal 
finite-dimensional filters, which in case of phase estimation are the well known phase track- 
ers, such as the phase locked loop (PLL), delay locked loop (DLL), angle tracking loops, and 
^ ' so on |14j . The MSEE in these phase trackers is asymptotically optimal [10], [T3] . 

The suboptimal phase trackers are known to lose lock (or slip cycles) |14| . The MTLL 
in these filters is simply the mean first passage time (MFPT) of the estimation error to the 
.£h ! boundary of the lock region. The MFPT from an attractor of a dynamical system driven 

^ ■ by small noise has been calculated by large deviations and singular perturbation methods 

[T5] . [15] . [T7j . and in particular, for the PLL [18]. The MTLL in particle filters for phase 
estimation was found in [19]. It has been found recently that minimizing the MNE leads 
to a finite, yet much longer MTLL than in the above mentioned phase estimators [20], [21] . 
This raises the question of designing a causal (or noncausal) phase estimator with maximal 
MTLL. 

The MTLL is the fundamental performance criterion in phase tracking and synchroniza- 
tion systems. Thus, for example, a phase tracking system is considered locked, as long as 
the estimation error e(t) = x(t) — x(t) is in (— 7T, 7t). When the error exceeds these limits, 
the estimation is said to be unlocked, and the system relocks on an erroneous equilibrium 
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point, with a deviation of 2ix. Another example is an automatic sight of a cannon. The sight 
is said to be locked on target if the positioning error is somewhere in between certain limits. 
Similar problems, in which the maximization of exit time is an optimality criterion, were 
considered by several authors [TJ, [8], [9]. In [7J, a simpler filtering problem is considered, in 
which the error e(t) is measured, rather than the state variable x(t). It is solved under the 
further assumption of a linear measurement inside a domain. In [8], [9] the state process is 
controlled through its drift, rendering it a control rather than a filtering problem. 

In this paper we show that for small noise the maximum MTLL filter is Bellman's MNE 
filter H|. It follows that the result of [ZE] for the MTLL of the optimal MNE phase filter, is 
asymptotically an upper bound for any other filtering scheme. In view of the results of |21j . 
the potential gain of the optimal MNE filter over the first order EKF-PLL is 12 dB. 

2 Formulation 

An important class of filtering problems with small measurements noise can be reduced to 
the model of a diffusion process 

dx(t) = m(x,t) dt + ea dw(t), (1) 

measured in a noisy channel 

dy{t) = h(x,t)dt + epdv(t), (2) 

where m(x,t) and h(x,t) are possibly nonlinear, continuous functions. The processes w(t) 
and v{t) are independent standard Brownian motions, and e is a small parameter. If m(x, t) 
is a linear function and the noise in ([T]) is not small, an appropriate scaling of time and 
dependent variables scales the small measurements noise into the diffusion equation as well, 
giving the canonic system (jTJ), ([2]) [20]. The optimal filtering problem is to find a causal 
estimator x(t) of x(t), given the measurements y\ = {y(s) : < s < £}, such that the mean 
first time the error signal, 

e(t) = x(t) -x(t), (3) 

leaves a given lock domain L C K, is maximal. More specifically, for any adapted function 
x(t) G C(1R + ) (measurable with respect to the filtration generated by y(t)), we define an 
error process by (J3J) and the first time to lose lock by 

r = inf {t : e(t) G dL} . (4) 

The optimal filtering problem is to maximize E[r \ y$] (see definition ffl6|) below) with respect 
to all adapted continuous functions x(t). For example, if h(x, t) = sinx in a phase estimation 
problem, then L = (— n, tt) and lock is lost when e(t) = ±7T. 

We can rewrite the model equations dlj, (j2J) in terms of the error process e(t) as 

de(t) = M £ (e(t),t)dt + sadw(t) (5) 

dy(t) = H £ (e(t),t)dt + epdv{t) (6) 
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where 

M £ (e(t),t) = m(x(t) + e{t)) - x(t) 

Hi{e{t),t) = h(x(t) + e(t)), 

and the filtering problem is to find x(t), such that E[r \ y^] is maximal. 

The survival probability of a trajectory (e(t),y(t)) of fl5]) with absorption at dL and 
can be expressed in terms of the pdf p(e, y, t £, 77, s) of the two-dimensional process with an 
absorbing boundary condition on dL. It is the solution of the Fokker-Planck equation (FPE) 

9p(e,y,t \£,7),s) dM £ (e,t)p(e,y,t\(,r], s) dH £ (e,t)p(e,y,t\ £,r),s) 



dt de dy 

e 2 a 2 d 2 p(e, y, 1 1 f , 77, s) e 2 p 2 <9 2 p(e,y, £ I £, 77, s) 



(7) 



2 de 2 2 dy 2 

for e,(£L, 7/, 77 G R, with the boundary and initial conditions 

p(e,y,t\£,r],s) = for e G <9L, 7/ G R, £ G L, 77 G R (8) 

p(e,j/,a|^77,s) = 5{e-^y-rj) for e G L, y G R, £ G L, 77 G R. (9) 

The pdf is actually the joint density and probability function p(e, y,t\£, 77, s) = Pr{e(£) = 
e, y{t) = y,r > £ | £, 77, s} and thus the survival probability is 

Pr{r > £| £,7j,s} = S e r. )>y (.)(t) = p(e,y,t\£,r],s) dedy, (10) 



and it decays in time. 



3 Simulation with particles 

To simulate the filtering problem on a finite interval < £ < T, we discretize ([I]), ([21 on a 
sequence of grids 

T 

U = iAt, i = 0,l,...,N, At 



N J ' 

and define discrete trajectories by the Euler scheme 

xn(U+i) = x N (ti) + Atm(x N (ti),ti) + ea Aw(ti) (11) 

= Vn{U) + Ath(x N (ti),ti) +epAu(t i ) J (12) 
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for i = 0, 1, . . . , N—l, where Aw(ti) and Av(ti) are independent zero mean Gaussian random 
variables with variance At. The discretized version of (151) is 

e N (t i+1 ) = e N {t i )+AtM & {e N {t i ),t i ) + eaAw{t i ) (13) 



VN{t 



y N (U) + AtH £ (e N (U),U) +epAu(U 



(14) 



Given an observed trajectory {yN(ti)}iL , we sam ple n trajectories {{£j,jv(£i)}i!=o}" = i> ac ~ 
cording to the scheme ffTTj) . which produce error trajectories {{e^jv^i)}^}™-!' anc ^ deter- 
mine their first exit times from L, denoted {t J; j\t}" =1 (we set tj jN = T if {e^jv^j)}^ does 
not exit L by time T) [22], [23], [21], [23], [26]. The conditional MTLL is defined on the 
ensemble by 



E 



r N A T | <{ y N (ti), i = 0, 1, . . . 

r jiA r A T 



At 



AiV 



(15) 



y^(Tj,jy A T) exp < 

3=1 



e 2 p 2 



At 

E 

fc=0 



H(e^ N (t k _i),t k -i)Ay KN - -if 2 (e J;JV (4-i), £ fe -i)At 



^exp 

3=1 



r jiA r A T 



At 

E 

fc=0 



^ (ej,jv(tfc-i),tfe-i)Ay fc) jv - -if 2 (e iiA r(t fc -i), t k -i)At 



We define 



^[ r l?/o] = lim lim um ^ ^AT| ^yjv(ti), i = 0, 1, 



I At. 



AN 



(16) 



The conditional MTLL i?[r | is a random variable on the cr-algebra J 7 ^ = ^J^t, where 



JF t is the cr-algebra generated by the measurements process y(-) up to time t Our purpose 
is to find x(t) that maximizes E[t | t/J] in the class of continuous adapted functions. 



4 The joint pdf of the discrete process 



The pdf of a trajectory of (ejv(t), VN{t)) is the Gaussian 



N 



PN[ei, e 2 , 



ejv; 3/1,3/2 



) • • • ; ifiv , 



£l, t 2 , • • • , £jv) — ] 



fc=l 



exp 



2e 2 At 



2ire 2 paAt 



, (17) 
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where the exponent is the quadratic form 

B k (x k , x k -x) = [x k - x k ^ - Ata k ^] T B [x k - x k -j - Ata k ^} , 

such that 



x k 



Vk 



M & (e k ,t k ) 
Hx(e k , t k ) 



B 



a- 2 ' 
p~ 2 



The Wiener path integral [27], [28], [29], [30], [31 

P(e, V, 1 1 £, r), s) = 



lim / dei / de 2 ■ ■ ■ / de^-i 



dyi / dy 2 --- dy N ^ x x 



N 

n 

fc=i 



exp 



7V-1 
Bk(Kk,X k -i) 

2e 2 At 



N~l 



2ne 2 paAt 



18) 



with 6 7v = e, = y, eo = £, yo = ? 7, is the solution of the FPE ([7]) with the boundary and 
initial conditions (18]) and (19]). 



(19) 



The pdf ( 1171) can be written as 

p N (ei, e 2 , • • • , e N ; yi, y 2 , • • • , 2/jv; *i, *2, • • • , *jv) 



JV 

n 

k=l 



1 



.A/27rAt£:o- 
1 



exp 



[e k - e fc _i - AtM 4 (e fc _ 1 , t 



fc-ij 



2£ 2 <X 2 At 

1 



exp 



/v exp 

n— 

^=1 



H £ (e k - 1 ,t k _ 1 )(y k -yk-i) - n ~ 2 2 fl|(e fc -i, 4-0 At 



2e 2 p 2 



(yfc - 2/fc-i) ; 

2e 2 p 2 At 



V2nAiep 



where, by the Feynman-Kac formula [27], [28], [29], [30], [31], the first product gives in the 
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limit the function 



exp 



cxp 



lim / dei / de 2 ■ • • 



f n 



N-l 



[e k - e fc _i - AtM a (e fc -i,t fc - 
2e 2 a 2 At 



x 



^H £ (e k ^ u t k _ 1 )(y k - y fc _i) - ^-^#?(e fc _i, £ fc _i)At 



which is the solution of the Zakai's equation in Stratonovich form [32 

d S (p(e,t,p) = 

y9(e,t)if £ (e,£) 



[ Mg(e, ^(e, t) ] e + ^VV(e, t ] ee - > dt + 



d s y{t), 



with the boundary conditions 

<p(e, £,/?) = for e e 9L. 

Therefore the joint density 

p N (e N ,t N ; yi,y 2 , ...,y N ) = 



(20) 
(21) 



Pr{e N (t N ) = e N ,T> t; y^h) = y u y N {t 2 ) = y 2 , ■ ■ ■ , 2/jv(*jv) = Vn} 
can be written at t — t N , = e as 



p N (e,t; 2/i, 3/2, ■ • .,Un) 



[(p{e, t, p) + o(l)] TT / -J_ exp 
^ V27rAtep 



2e 2 p 2 At 



where o(l) — * as iV — >• oo. Equivalently, 

Piv(e,t; 2/i,2/2, • • • , Vn) 



<p(e,t,p) 



N 



^v/^Afc/^l 2eVAt 



(22) 



(23) 



which can be interpreted as follows: tp(e, t, p) is the joint conditional density of ejv(i) and r > 
t, given the entire trajectory {yN(U)}iLo, however, the probability density of the trajectories 
{Vn(U)}1 , 

(y k - ?/fc-i) 2 ^ 



A' 



Pn 



(vl) = I 



k=l 



exp 



2e 2 p 2 At 
V2nAtep 
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is Brownian, rather than the a priori density imposed by (151), ([6]). 

Now, 

Pr{r > t N , y N (ti) = yi,y N (t 2 ) = 2/2, • • • , Vn^n) = Vn} = 



Pr{r > t N | 2/jv(*i) = 2/1, Iforfcj) = 2/2, ■ ■ ■ , 2/jv(*jv)} x 

Pr{2/iv(*i) = 2/i, I/j\r(*2) = 2/2, ••• , = 2/tv}, 

which we abbreviate to 

Pr{r > t, yl) = Pr{r > t 1 yl}p N {vl), (24) 

where the density pjv(yo) = Pr {2/Ar(ti) = 2/i, 2/w(^) = 2/2, • • • , 2/at(£/v) = Vn} is defined by the 
system ( fill) . ( |T2i) . independently of 

We now use the abbreviated notation (|24|) to write 

or ^ ., t1 Pr{ r > t, 2/jv(*i) = UuVn(1*) = 2/2, ••• ,2/jv(*jv) = 2/tv} 

Prjr > t 2/o) = r~F\ 

PM2/o 



p N (e,t;y 1 ,y 2 , . . . , 2/jy) 



= ^TTT / Me,*,p) + o(l)}de. (25) 
As N ^ oo, both sides of eq. (l25j) converge to a finite limit, which we write as 

Pi{T>t\yl} = a(t) J <p(e,t)de, (26) 

where 

a(t) = hm 

is a function independent of x(t). 

Next, we show that E[r \ z/J], as defined in (fT5l) . ffl6|) . is given by 

roo 

E\r\yl}= / Pi{r>t\y t }dt. (27) 
./o 

Indeed, since Pr{r > t \ y$} — > exponentially fast as t — > oo, we can write 
/ Pr{r>t|^}dt= lim / tdPr{r < 1 1 z/q} 

and 

/ tdPr{r <t\y* } = lim V zAtA Pr{r < zAt | y l At }, 

J0 8=1 
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where 



A Pr{r < iAt \ y* At } = Pr{r < iAt \ yf 1 } - Pr{r < (i - l)At \ y^~ L)ax }. 



(i-l)At-| 



Now, we renumber the sampled trajectories ej^iU) i n the numerator in (|T5l) according to 
increasing r^jy, so that in the new enumeration r^jy = zAi. Then we group together the 
terms in the sum that have the same r, m and denote their sums m,- at, so that (1 151) becomes 



E 



t n A T | <^ j/7v(ti), i = 0, 1, . . . 



tat 
At 



A AT 



i=l 



AT 

8=1 



Finally, we identify 



APr{r < zAi | 



mj,N 

N 
i=l 



where o(l) ^ as AT — > oo. Hence (|27[) follows. Finally, we identify 



A Pr{r < i At \ y l Q m } 



N 



i=l 



where o(l) — > as A" — > oo. Hence (|27[) follows. 



(2? 



4.1 Asymptotic solution of Zakai's equation and the optimal filter 

For small e the solution of (1201 with the boundary conditions (12 ip is constructed by the 
method of matched asymptotics [33], [16], [17]. The outer solution is given by large deviations 
theory [T3], [M], [M], [33] as 



V outer (e,i) = exp 



where 

^(e,t) 

and 



inf 

e(-) 6 <*([(), *])■/<> 



e(s)-M 4 (e(a),a) 



y(s) - H & (e(s),s) 



ds, (29) 



£([0,*]) = {e(-)eC\[0,t}) :e(0) = e} 
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We denote by e(t) the minimizer of the integral on the right hand side of eq. (l29l) . The outer 
solution y 9 outcr( e ) ^) does not satisfy the boundary conditions (12TI) . so a boundary layer 
correction k(e,t,e) is needed to obtain a uniform asymptotic approximation, 



tp(e, t) ~ Vuniform(e, t) = ^ uter(e, t, p)k(e, t, e) = exp 



V(e,t) 



k(e, t, e). 



The boundary layer function has to satisfy the boundary and matching conditions 
k(e,t,e) = for e G <9L, lim k(e, t, e) = 1 for e G L, 

e— >0 

uniformly on compact subsets of the interior of L. 
Since the survival probability is 



Pr{r>i|y$} = ^a(t)exp j- 



the MTLL, according to ( 1271) . is given by 



a(t) exp 



.^2^1 J, k(e,t,£) de, 



(30) 



(31) 



(32) 



In view of ([3]), the minimizer e(t) of the integral on the right hand side of (1291) can be 
represented as e(t) = x(t) — x(t), where x(t) is the minimizer of the integral 



= inf 

x{.)eCl{[0,t])Jo 



tAf 



x{s) — m{x{s), s) 



a 



y(s) - h(x(s),s)' 
P 



ds, (33) 



where r = inf{£ : x(t) — x(t) G dL) and 

C 1 x ([0 1 t]) = {x(-)eC 1 ([0,t]) :x(0) = x}. 
Writing ij)(e,t) = ty(x,t) and k(e,t,e) = K(x,t,e), we rewrite ( l32l) as 



E[r | ?/ T ] 



a(t) exp <| ^ — [• £, e) dx dt. 

•/£+&(*) I e 



(34) 



The integral in (1341) is evaluated for small e by the Laplace method, in which the integrand 
is approximated by a Gaussian density with mean x(t) and variance proportional to e 2 . 
It is obviously maximized over the functions x(t) by choosing x(t) so that the domain of 
integration covers as much as possible of the area under the Gaussian bell. If L is an interval, 
then the choice x(t) = x(t) is optimal. We conclude that for small noise, the minimum noise 
energy filter x(t) is asymptotically the maximum MTLL filter. 
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5 Discussion 



The main result of this paper is a proof that for small noise, the minimum noise energy filter 
maximizes the mean time the estimation error stays within a given region, e.g., maximizes 
the mean time to lose lock in problems of phase tracking and synchronization. The MNE 
filter is not finite-dimensional, however finite discrete approximations, such as Viterbi-type 
algorithms [37], J3S], can give arbitrary accuracy. The practical aspects of finding the true 
MNE filter, or otherwise adequate approximations for it, was partially dealt with in [2Tj and 
still remains an interesting issue for further studies. 

Katzur et. al. [36], and subsequently Picard [11] [12], have shown that for nonlinear, 
but monotone measurement functions, the MNE filter is to leading order identical to the 
extended Kalman filter. However, for measurement functions which are non-monotone, this 
is apparently not the case. Ezri [2U] and Fischler [21] have considered the problem of phase 
filtering and smoothing respectively, in which the stochastic phase process x(t) is measured 
in a low noise channel by the vector function h(x) = [sin(x), cos(x) ] T . They show that 
there is a huge gap between the MTLLs of the extended Kalman filter (smoother) or particle 
filter, and the MNE filter (smoother), respectively. 

The great advantage of the MNE filter in the case of phase estimation is explained by 
the observation that finite-dimensional approximations to the MAP or minimal MSEE filters 
(the EKF or the finite dimensional filters of Katzur [36]), do not capture large deviations of 
the signal or of the measurements noise. They are optimal only near local maxima of the a 
posteriori probability density. The MNE filter, in contrast, is a global MAP estimator and 
can track large deviations. Thus, it is less vulnerable to loss of lock phenomena, relative to 
the above mentioned filters. 

Acknowledgment: The authors thank B.Z. Bobrovsky, O. Zeitouni, D. Ezri, B. Nadler 
and A. Taflia for useful discussions. 
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